file1=read.table("ASM_Bayes_meta_anno_m1.txt",head=T,sep="\t")
file2=read.table("E:/1.甲基化分析/ASM/ASM_snp-onlyWGS/ASM_Bayes_Meta_Chisq_anno_num2_m1_new.txt",head=T,sep="\t")
file22=data.frame(file2[,1:14],file2[,56:63])
file11=data.frame(file1[,1:14],file1[,28:35])
header=c("unitID",paste(names(file11)[2:22],"5hmc",sep="_"))
names(file11)=header
header=c("unitID",paste(names(file22)[2:22],"5mc",sep="_"))
names(file22)=header
file3=merge(file11,file22,by="unitID")
###用beta0来比较
c=data.frame(unitID=file3$unitID,file3[,19:22],file3[,40:43],snpid=file3$avsnp150_5hmc)
c=c[order(c$DC_case_mean_5hmc),]
a=data.frame(c[,1:5],snpid=c$snpid)
b=data.frame(unitID=c$unitID,c[,6:10])
b$num=1:17
b$group="5mc"
a$num=1:17
a$group="5hmc"
header=c("unitID","DC_case_pvalue","DC_case_mean","DC_case_down","DC_case_up","snpid","num","group")
names(a)=header
names(b)=header
cc=rbind(a,b)

ggplot(cc,aes(y=cc$DC_case_mean,x=cc$num))+scale_x_continuous(breaks = cc$num,labels = cc$snpid)+
geom_errorbar(aes(ymin=cc$DC_case_down,ymax=cc$DC_case_up,color=group))+geom_point(aes(color=group),size=1)+
labs(x="snpID",y="β0")+geom_hline(yintercept = 0,color="green",size=1)+coord_flip()+theme_bw()

###用case的var allele reads之和来比较
id=as.character(file3$unitID)
row.names(file1)=file1$unitID
row.names(file2)=file2$unitID
file11=file1[id,]
file22=file2[id,]

sel1=c("X2B_X1T","M8_M7","M6_M5","M2_M1","M48_M47","M50_M49")
sel2=c("M30_M29","M26_M25","M35_M36","M28_M27")
sel1=unlist(strsplit(sel1,"_"))[seq(2,2*length(sel1),2)]
sel2=unlist(strsplit(sel2,"_"))
sel=c(paste(sel1,"_reads2",sep=""),paste(sel2,"_reads2",sep=""))
file11=file11[,c("unitID","avsnp150",sel)]
file22=file22[,c("unitID","avsnp138",sel)]
file22[file22=="."]=NA
file12=file11[,3:16]
file23=file22[,3:16]
###计算sum
file11$sum=0
for (i in 1:17){
#sels=which(is.na(file12[i,]))
file11[i,]$sum=sum(file12[i,],na.rm=T)
}

file22$sum=0
for (i in 1:17){
file22[i,]$sum=sum(file23[i,],na.rm=T)
}

a=data.frame(unitID=file11$unitID,snpid=file11$avsnp150,sum=file11$sum)
a=a[order(a$sum),]
a$group="5hmc"
a$num=1:17
c=data.frame(unitID=a$unitID,num=a$num)

b=data.frame(unitID=file22$unitID,snpid=file22$avsnp138,sum=file22$sum)
b$group="5mc"
b=merge(b,c,by="unitID")
cc=rbind(a,b)

ggplot(cc,aes(x=cc$num,y=cc$sum))+scale_x_continuous(breaks = cc$num,labels = cc$snpid)+
geom_point(aes(color=group),size=1)+geom_line(aes(color=group),size=1)+theme_bw()+theme(axis.text.x = element_text(angle = 60,hjust=.9))
###均值
file11$mean=0
file22$mean=0
for (i in 1:17){
#sels=which(is.na(file12[i,]))
file11[i,]$mean=mean(unlist(file12[i,]),na.rm=T)
file22[i,]$mean=mean(unlist(file23[i,]),na.rm=T)
}

a=data.frame(unitID=file11$unitID,snpid=file11$avsnp150,mean=file11$mean)
a=a[order(a$mean),]
a$group="5hmc"
a$num=1:17
c=data.frame(unitID=a$unitID,num=a$num)

b=data.frame(unitID=file22$unitID,snpid=file22$avsnp138,mean=file22$mean)
b$group="5mc"
b=merge(b,c,by="unitID")
cc=rbind(a,b)
ggplot(cc,aes(x=cc$num,y=cc$mean))+scale_x_continuous(breaks = cc$num,labels = cc$snpid)+labs(x="snpID",y="mean")+
geom_point(aes(color=group),size=1)+geom_line(aes(color=group),size=1)+theme_bw()+theme(axis.text.x = element_text(angle = 60,hjust=.9))

###VAF
sel1=c("X2B_X1T","M8_M7","M6_M5","M2_M1","M48_M47","M50_M49")
sel2=c("M30_M29","M26_M25","M35_M36","M28_M27")
sel1=unlist(strsplit(sel1,"_"))[seq(2,2*length(sel1),2)]
sel2=unlist(strsplit(sel2,"_"))
sel=c(paste(sel1,"_var_freq",sep=""),paste(sel2,"_var_freq",sep=""))
file11=file11[,c("unitID","avsnp150",sel)]
file22=file22[,c("unitID","avsnp138",sel)]
file22[file22=="."]=NA
file12=file11[,3:16]
file23=file22[,3:16]

t12=t(file12)

file11$mean=0			#算得排序的序号
file22$mean=0
for (i in 1:17){
#sels=which(is.na(file12[i,]))
file11[i,]$mean=mean(unlist(file12[i,]),na.rm=T)
file22[i,]$mean=mean(unlist(file23[i,]),na.rm=T)
}

a=data.frame(unitID=file11$unitID,snpid=file11$avsnp150,mean=file11$mean)
a=a[order(a$mean),]
a$num=1:17
c=data.frame(unitID=a$unitID,num=a$num)

name=row.names(file12)


for (i in 2:ncol(t12)){
part=data.frame(t12[,i])
names(part)="VAF"
part$unitID=name[i]
part=part[!is.na(part[,1]),]
part$group="5hmc"
part=merge(part,c,by="unitID")
final=rbind(final,part)
}

t23=t(file23)
name=file22$unitID

for (i in 2:ncol(t23)){
part=data.frame(t23[,i])
names(part)="VAF"
part$unitID=name[i]
part=part[!is.na(part[,1]),]
part$group="5mc"
part=merge(part,c,by="unitID")
finals=rbind(finals,part)
}

ggplot(finale,aes(x=finale$VAF,y=finale$num))+scale_y_continuous(breaks = finale$num,labels = finale$unitID)+labs(x="VAF",y="snpID")+
geom_point(aes(color=group),size=3)+theme_bw()

